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Global Gyrokinetic Simulation of Tokamak Edge Pedestal Instabilities 
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Global electromagnetic gyrokinetic simulations show the existence of near threshold conditions 
for both a high-n kinetic ballooning mode (KBM) and an intermediate-n kinetic version of peeling- 
ballooning mode (KPBM) in the edge pedestal of two DIII-D H-mode discharges. When the magnetic 
shear is reduced in a narrow region of steep pressure gradient, the KPBM is significantly stabilized, 
while the KBM is weakly destabilized and hence becomes the most-unstable mode. Collisions 
decrease the KBM's critical /3 and increase the growth rate. 

PACS numbers: 52.25.Fi, 52.35.-g, 52.55.Tn, 52.65.-y 



Present-day tokamak fusion experiments achieve high 
performance with a narrow edge particle and energy 
transport barrier at the plasma boundary, called an 
"edge pedestal." In the edge pedestal, the plasma den- 
sity and temperature gradients build up and eventually 
drive electromagnetic instabilities at intermediate and 
high toroidal mode numbers (n). Understanding the un- 
derlying physics of pedestal instabilities from first prin- 
ciples is critical for predicting the performance of future 
large experiments such as ITER [1[. The pedestal can 
be crudely characterized by the height and width of the 
pressure profile. Snyder et al. [2|, |3[ developed a model 
that successfully predicts the experimentally observed 
height and width of the pedestal by combining the lin- 
ear threshold of two electromagnetic plasma instabilities 
in the pedestal region: the intermediate-n magnctohy- 
drodynamic (MHD) "peeling-ballooning" mode (PBM) 
and the high-n kinetic ballooning mode (KBM). While 
the PBM has been well parameterized using ideal MHD 
theory, it has never, until now, been verified using more 
realistic kinetic calculations. The KBM, on the other 
hand, requires a kinetic model. 

The KBM threshold is known to be important in core 
turbulence simulations Ij,l5|- Transport levels in core cal- 
culations become very high when approaching the KBM 
threshold. The KBM threshold has not been cle arly 
identified in previous gyrokinetic edge simulations [6J-llO| . 
However, experimental evidence of KBMs exist in DIII-D 
H-mode (High confinement mode) and quiescent H-mode 
experiments [ll|, [l^, and the observed profiles closely 
correspond to a simplified calculation of KBM critical- 
ity over a wide range of par ameters [3|, ll3|. Recently, 
GS2 simulations of MAST IJ, ll5| indicated an even par- 
ity mode in the steep gradient region of the pedestal of 
the spherical tokamak and identified the mode as a KBM. 
The main signature of a KBM - an electromagnetic mode 
with a /3 threshold, has still yet to be demonstrated, until 
this Letter. 

We present a self-consistent picture of both the KBM 
and PBM in the pedestal, and describe the conditions 



when each of them dominates. We show that the H-mode 
pedestal, just prior to the onset of observed Edge Local- 
ized Mode (ELM) instabilities, is very near the KBM 
threshold in global gyrokinetic simulations. In addition 
to the high-71 KBM, an intermediate-n electromagnetic 
mode is unstable and we identify it as a kinetic ver- 
sion of the MHD peeling-ballooning mode. Using the 
gyrokinetic 5f particle- in-cell code GEM [l6|, ll7| with 
electron-ion collisions, we study the global linear sta- 
bility of H-mode pedestal profiles from two DIII-D ex- 
periments: discharge 136051 that has been previously 



reported [11[ with characteristics of KBM, and another 
discharge 132016. Calculations using these "original" 
profiles show two types of instabilities: an intermediate- 
n mode that propagates in the electron diamagnetic di- 
rection in the plasma frame (we will call this mode the 
"kinetic peeling ballooning mode", KPBM) and a high- 
n. low frequency mode that mostly propagates in the 
ion direction (we refer to this mode as the "ion mode"). 
These two modes are driven by the pressure gradient. 
While the ion mode has a ballooning structure, it lacks 
an important property of the KBM: the /3 parameter 
scan should show a sudden change of real frequency cor- 
responding to strong increase of growth rate (see, e.g. 



Refs [18|, [I9j ) . Additionally, this ion mode is subdomi- 



nant to the KPBM. The KPBM is very sensitive to the 
shape of the q profile and can only be seen in global sim- 
ulations and not in flux tube simulations. If we slightly 
manipulate the magnetic shear near the steep pressure 
gradient region of the pedestal, i.e., very locally flatten 
the q profile, the KPBM would be significantly stabilized 
and the ion mode would be weakly destabilized and begin 
to show clear KBM characteristics. These results indi- 
cate that an improved pedestal model should include, in 
detail, any corrections to the bootstrap current |20l422| . 
and any other equilibrium effects that might reduce the 
local magnetic shear 23|, l2J] . 



The two experimental profiles are shown in Fig. [T] The 
simulation box covers the 0.899 < pn < 0.999 region in- 
side the separatrix, where pN is the normalized radius. 
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FIG. 1: (Color online). The experimental profiles of (a) Ti, 
(b) Te, (c) density n and (d) original and flat safety factor q 
for 136051 (solid lines) and 132016 (dash-dotted lines). 



Fixed boundary conditions are applied and the density 
and temperature profiles are smoothed at the boundaries 
in simulation. The magnetic equilibria are parameterized 
using Miller equilibrium [71, |25| . The magnetic equilib- 
rium used for 136051 did not include corrections for the 
bootstrap current. The equilibrium for 132016 included 
corrections using the Sauter model |2l| . The simulation 
domain grid is 64 x 32 x 32, with 64 cells along the radial 
direction. The time step is Ai = 1/0.^ where fli is the 
proton gyrofrequency calculated at top of the pedestal. 
There are 1048576 particles per species with realistic deu- 
terium to electron mass ratio. Figure [ijd) shows the 
flattened q profiles as well. 

Figure [2] scans instabilities with mode number of 7 < 
n < 70 for the two profiles with the original q profiles 
from experiments. The corresponding kyPD at the cen- 
ter of the simulation box is in the range of [0.102, 1.02] 
for 136051 and [0.109, 1.09] for 132016, where po is deu- 
terium gyroradius. The two discharges exhibit quite sim- 
ilar trends. From Fig. ^a.) there are clearly two types of 
instabilities: intermediate-n {n < 21) modes and high-n 
modes. Both modes appear to propagate in the elec- 
tron diamagnetic direction here indicated by their posi- 
tive real frequencies. However, there is a Doppler shift 
caused by the radial electric field E^ here. In simula- 
tions without Er shown later in Figs. [5] and [6l the high- 
n instability propagates in the ion direction while the 
intermediate-n instability still in the electron direction. 
The experiment of 136051 [111 has found two bands of 
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FIG. 2: (Color online). The results of (a) linear real frequency 
Ur and (b) growth rate 7 for the two discharges with original 
q profiles. The effect of collisions is also shown. 



density fluctuations, with an ion band at 50 - 150 kHz 
and an electron band at 200 - 400 kHz. Here for 136051, 
in the "laboratory" frame, the high-n instability has a 
frequency equivalent to 160 kHz, quite close to the fre- 
quency of the ion band found in experiment. However, 
unlike in the experiment, this mode is not the dominant 
instability here. As shown in Fig.[2jb), the intermediate- 
n electron instability has a much higher growth rate. Its 
frequency is around 800 kHz here, about twice that of 
the electron band of the experiment. 

The electrostatic potential (j){x, y) of these two insta- 
bilities are shown in Fig. |3|a,b) for 136051, with x and y 
corresponding to the radial and toroidal direction in the 
field-line-following coordinate, respectively. Figure [3l[c) 
also shows the temperature and density gradients, repre- 
sented by R/Ln and R/ Lt, where R is the major radius 
and L~^ = d\nn/dr and L^^ = d\nT/dr. Note that the 
gradients are zero at the boundaries because the profiles 
are smoothed. The x axes of the contour plots corre- 
sponds to the radius of Fig.[31Jc). Both instabilities peak 
in the steep gradient region, indicating they are driven by 
pressure gradients. Both instabilities also have a largely 
even parity structure. However, the intermediate-n mode 
appears to have a "tail" tilted towards the top of the 
pedestal; the dominant structure of the high-n mode has 
a tail tilted towards the separatrix. In fiux tube simula- 
tions we find ITG is the dominant instability on top of 
the pedestal, but its growth rate is weaker than that of 
the high-n mode in the steep gradient region and there- 
fore ITG is never dominant in global simulations. The 
possibility of trapped electron mode [26[ is excluded be- 
cause collisions don't decrease the linear growth rate. 




FIG. 3: (Color online). The mode structure (electrostatic 
potential contour plots) of the intermediate-n (a) and high-n 
(b) modes of 136051, with the original q profile, and compared 
to the temperature and density gradients (c). 



The results of the high-n instability (ion mode) here 
agree with our flux tube simulations. In simulations 
with another profile 131997 [8[, a similar mode was 
found and the results agree with electrostatic simula- 
tions of GTC [10| and the flux tube eigenmode results 
of GYRO [9|, although the mode has a positive real fre- 
quency for kypo < 0.35 without Er and was thus iden- 
tified as an "electron mode" . The intcrmediate-n mode, 
which we now refer to as KPBM, is not observed in flux- 
tube simulations and is very sensitive to the q profllc. 
Measuring the pedestal q profllc is a difficult experimen- 
tal challenge. Thus, experimental values for the q-profile 
are usually obtained from application of bootstrap cur- 
rent models to measured pedestal density and tempera- 
ture profiles. In fact, the bootstrap current can vary from 
the generally used Sauter model [22] , and there is signif- 
icant uncertainty in the measured gradients required to 
calculate the bootstrap current. In previous simulations 
of discharge 98889 [27| , which has a near-zero magnetic 
shear in a region across the steep gradient area [2J], the 
KPBM is not present. We now fiattcn the q profiles in 
a very small region in the two discharges as shown in 
Fig. [IJd). In doing so, we run flux tube simulations flrst 
to flnd a position in the steep gradient region that is lo- 
cally most unstable, and then change the q profile at that 
location with a zero magnetic shear. 

Figure [4] shows the results with the "flat" q profiles. 
The intermediate- ri KPBM is significantly stabilized and 
the high-n ion mode, which we now identify as KBM, 
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FIG. 4: (Color online). Same as Fig. [2] but with flattened g's. 



now dominates. The flat q reduces the real frequency for 
ri = 14 and n = 21, making the frequencies comparable 
with KBM. Collisionality further suppresses the high fre- 
quency of the n = 7 mode and reduces its growth rate. In 
addition, the (t>{x, y) mode structures of the modes (not 
shown) arc also changed from Fig. [3l the tilted structure 
is reduced and the modes exhibit the more typical even 
structure. If we reduce the magnetic shear gradually in- 
stead of using a flat q here, the growth rate of the KPBM 
is also reduced gradually. 
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FIG. 5: (Color online). Simulations of 136051 at a high /3 and 
without Er. Results with the flat q profile are compared to 
that with the original q, including the effect of collisionality. 



We can see both the effects of the flat q and collision- 
ality more clearly in runs at higher values of /3, where 
the modes are more unstable. Figure [5] shows results of 
136051 with twice the experimental /3 for both the orig- 



inal and flat q profile. Er is removed to eliminate the 
Doppler shift. Although E^ is generally believed to be a 
stabilizing factor at pedestal, here in linear simulations 
it is found to be destabilizing for the KPBM, thereby the 
growth rates in Fig. [SJa) are smaller than in Fig. [IJb). 
The KBM now has a negative real frequency. It becomes 
obvious that the flat q significantly stabilizes the KPBM 
and reduces its real frequency, leaving the KBM domi- 
nant. The KBM is moderately destabilized by the flat 
g, with its real frequency almost unchanged. Collisions 
reduce and even suppress the frequency of the KPBM, 
and are slightly destabilizing for the KBM. For n > 70 
the growth rate may rise again, but we restrict this study 
to kyPi < 2 where our gyrokinetic simulations arc valid. 

Figure [6] shows a /3 scan for the KBM for both ex- 
perimental profiles. In simulations with both the slightly 
flattened q and collisions, the instabilities of both profllcs 
display the standard KBM features. As f3 increases, the 
growth rate remains low (electrostatic), and then after 
passing the critical f3, the growth rate strongly increases, 
with a corresponding sudden change in the real frequency 
with phase velocity in the ion direction. Since the exper- 
imental /3 (/3x = 1) is well above the critical /3, the KBM 
is unstable, and is indeed the dominant instability, which 
in turn would limit the pressure gradients of the pedestal. 
The effect of collisions is to make the critical f3 smaller. 
In simulations with the original q profiles, there's not the 
characteristic sudden change in the real frequency, but 
the mode could still be identified as a modified KBM 
based on the linear results presented here. 
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Resistive MHD studies of the PBM mode by Zhu et al. 
J23| have shown that a flat q stabilizes the PBM, the 
same property we see here with the KPBM. Besides be- 
ing sensitive to the oprofile, the PBM is also electromag- 
netic, and in Ref. |23| it has n < 11, a similar range of 
mode number of our KPBM. The unique property of the 
KPBM here is that it has a phase velocity in the electron 
diamagnetic direction while resistive MHD would show a 
near zero real frequency. 

Theoretically, the effects of magnetic shear (charac- 
terized by s = d\Tiq{p)/d\'np) and pressure gradient 
(characterized by a = q^R^/Lp, where Lp = p/Vp is 
the pressure scale length) on the stability of ballooning 
modes have been studied extensively. For core plasmas 
it is shown that the ideal MHD ballooning modes could 
be destabilized by a near-zero magnetic shear, as in the 
internal transport barrier [28[; while KBMs are found 
to be unstable with a negative [29[ or near-zero mag- 
netic shear [SO] ■ In the steep gradient region of the edge 
pedestal, however, the pressure gradient is so high and a 
is usually bigger than the stable threshold of typical core 
plasma s - a diagrams and the cases we studied here 
should be near the second stability region. The s - a 
analysis of edge plasmas could therefore be very different. 
Furthermore, our kinetic simulations suggest the gradi- 
ents of density and electron and ion temperatures have 
different destabilizing effects, which apparently cannot 
be represented by a single HMD parameter a. 

Transport codes [SlJ, MHD calculations 1231 and ex- 
perimental measurements 



13 



have all shown that the 
bootstrap current can flatten the q proflle. The kinetic 
linear stability of the edge pedestal is thus a subtle com- 
petition between the PBM and KBM as seen both here 
and previously in the EPED model [3,l3|. Both are driven 
by the pressure gradients and therefore limit the pedestal 
shape. If the magnetic shear is high, PBM is much more 
unstable than KBM. Reducing the magnetic shear stabi- 
lizes PBM and KBM becomes the dominant instability. 
It is thus important to incorporate accurate representa- 
tions of the bootstrap current and edge geometry, as both 
strongly impact magnetic shear. Additionally, better ex- 
perimental characterization of the edge q profile would 
help test and improve predictive models. 
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FIG. 6: (Color online). The /3 scan of 136051 at n = 49 and 
132016 at n = 56, with flat vs. original q profiles. "/3x" 
means the factor that is used to multiply the experimental /3. 
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